---
title: "Fixed Imaging Quantification"
output:
  html_document:
   code_folding: hide
---

```{r toolbox/importation, include=FALSE, warning = FALSE}
# Load packages
library(ggplot2)
library(moments)
library(dplyr)
library(emmeans)
library(ggsignif)
library(MASS)
library(lme4)
library(lmtest)
library(RColorBrewer)

# Import Data (csv from excel)
ens0 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Immunostaining Image Quantification/Fixed ENS Image Quant Combined Data v5.csv", header=TRUE, stringsAsFactors=FALSE)
cort0 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Immunostaining Image Quantification/Fixed cortex Image Quant Combined Data v2.csv", header = TRUE, stringsAsFactors = FALSE)
ens_tubb0 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Immunostaining Image Quantification/ens tubb3 v3.csv", header=TRUE, stringsAsFactors=FALSE)
cort_tubb0 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Immunostaining Image Quantification/cortex tubb3 v1.csv", header = TRUE, stringsAsFactors = FALSE)

# Define variables for analysis (e.g., group factor, dosage level)
ens0$Groups.f <- factor(ens0$group, levels = c("Control","But", "LPS", "But+LPS",
                                             "PFF", "PFF+But", "PFF+LPS", "PFF+But+LPS",
                                             "2XPFF", "4XPFF", "4XPFF+But", "6XPFF"))
ens0$attempt.no <- factor(ens0$attempt)
ens0$pff.qual <- ifelse(ens0$attempt.no == 6 & 
                            (ens0$group == "6XPFF" | 
                             ens0$group == "6XPFF+LPS" | 
                             ens0$group == "6XPFF+But" | 
                             ens0$group == "6XPFF+But+LPS"), 
                       1,0)
ens0$but.qual <- ifelse(ens0$attempt.no == 19 & (ens0$Groups.f == "But" | 
                                                  ens0$Groups.f == "But+LPS" | 
                                                  ens0$Groups.f == "PFF+But" | 
                                                  ens0$Groups.f == "4XPFF+But"), 1,0)
ens0$w.id <- factor(paste(ens0$attempt.no,"-",ens0$Groups.f,"-",ens0$well))
cort0$Groups.f <- factor(cort0$group, levels = c("Control","PFF","4XPFF"))
cort0$attempt.no <- factor(cort0$attempt)
ens_tubb0$attempt.no <- factor(ens_tubb0$attempt)
ens_tubb0$Groups.f <- factor(ens_tubb0$group, levels = c("Control","But", "LPS", "But+LPS",
                                             "PFF", "PFF+But", "PFF+LPS", "PFF+But+LPS",
                                             "2XPFF", "4XPFF", "4XPFF+But", "6XPFF"))
ens_tubb0$pff.qual <- ifelse(ens_tubb0$attempt.no == 6 & 
                            (ens_tubb0$group == "6XPFF" | 
                             ens_tubb0$group == "6XPFF+LPS" | 
                             ens_tubb0$group == "6XPFF+But" | 
                             ens_tubb0$group == "6XPFF+But+LPS"), 
                       1,0)
cort_tubb0$Groups.f <- factor(cort_tubb0$group, levels = c("Control","PFF","4XPFF"))
cort_tubb0$attempt.no <- factor(cort_tubb0$rat)
```

```{r unify dapi count cortex,  warning = FALSE}
# Import Nuclei count data (csv from excel)
cortdap.10 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Immunostaining Image Quantification/misc/DAPI Counting Analysis/Cortex Quant (except 14)/MyExpt_Nuclei.edit.csv", header = TRUE, stringsAsFactors = FALSE)
cortdap.20 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Immunostaining Image Quantification/misc/DAPI Counting Analysis/Cortex (Att 14)/MyExpt_Nuclei.csv", header = TRUE, stringsAsFactors = FALSE)

# Extract specific columns
cortdap.11 <- cbind(cortdap.10[,1],cortdap.10[,2],cortdap.10[,7],cortdap.10[,10],cortdap.10[,12],cortdap.10[,9])
cortdap.21 <- cbind(cortdap.20[,1],cortdap.20[,2],cortdap.20[,8],cortdap.20[,11],cortdap.20[,13],cortdap.20[,10])

# Merge two matrices
cortdap <- rbind(cortdap.11,cortdap.21)
cortdap <- data.frame(cortdap)

# Attach column names
colnames(cortdap) <- c("image.no",
                       "nuclei.no",
                       "group",
                       "rat.no",
                       "well",
                       "im.no")

# Create ID Variable
cortdap$id.maker <- paste(cortdap$rat.no,"-",
                          cortdap$group,"-",
                          cortdap$well,"-",
                          cortdap$im.no)

# Count number of nuclei (convert rows to individual nuclei to image)
cortdap.short <- data.frame(aggregate(data = cortdap,
                          nuclei.no ~ id.maker,
                          function(nuclei.no) length(unique(nuclei.no))))
cort0$id.maker <- paste(cort0$attempt,"-",
                          cort0$group,"-",
                          cort0$well,"-",
                          cort0$image.no)

cort.m2 <- merge(cort0,cortdap.short, all = TRUE)
cort <- merge(cort0,cortdap.short,by="id.maker")
```

```{r unify dapi count ens,  warning = FALSE}
# Import Nuclei count data (csv from excel)
ensdap.10 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Immunostaining Image Quantification/misc/DAPI Counting Analysis/ENS Quantification (except 14)/MyExpt_Nuclei.edit.csv", header = TRUE, stringsAsFactors = FALSE)
ensdap.20 <- read.csv("/Users/adambindas/Dropbox/Parkinson's disease ENS Paper/Data (Raw and Organized)/Immunostaining Image Quantification/misc/DAPI Counting Analysis/ENS (Att 14)/MyExpt_Nuclei.csv", header = TRUE, stringsAsFactors = FALSE)

# Extract specific columns
ensdap.11 <- cbind(ensdap.10[,1],ensdap.10[,2],ensdap.10[,7],ensdap.10[,10],ensdap.10[,12],ensdap.10[,9])
ensdap.21 <- cbind(ensdap.20[,1],ensdap.20[,2],ensdap.20[,8],ensdap.20[,11],ensdap.20[,13],ensdap.20[,10])

# Merge two matrices
ensdap <- data.frame(rbind(ensdap.11,ensdap.21))

# Attach column names
colnames(ensdap) <- c("image.no",
                       "nuclei.no",
                       "group",
                       "rat.no",
                       "well",
                       "im.no")

# Create ID Variable
ensdap$id.maker <- paste(ensdap$rat.no,"-",
                          ensdap$group,"-",
                          ensdap$well,"-",
                          ensdap$im.no)

# Count number of nuclei (convert rows to individual nuclei to image)
ensdap.short <- data.frame(aggregate(data = ensdap,
                          nuclei.no ~ id.maker,
                          function(nuclei.no) length(unique(nuclei.no))))

# Create ID Variable
ens0$id.maker <- paste(ens0$attempt,"-",
                          ens0$group,"-",
                          ens0$well,"-",
                          ens0$image.no)
ens_tubb0$id.maker <- paste(ens0$attempt,"-",
                          ens0$group,"-",
                          ens0$well,"-",
                          ens0$image.no)

ens <- merge(ens0,ensdap.short,by="id.maker")
ens_tubb.merge <- merge(ens_tubb0,ensdap.short,by="id.maker")
```

```{r Initial Inspection/Organization,  warning = FALSE}
#Subset datasets
ens2 <- subset(ens0, 
                 quality.removal == 0 & sat.qual.test == 0 & 
                 pff.qual == 0 & but.qual == 0 &
                 attempt != 4 & attempt != 14 & attempt != 6)
ens3 <- subset(ens2,
                 X6xpff == 0 & X6xpff.but == 0 & X6xpff.lps == 0 & X6xpff.but.lps == 0 & 
                 X2xpff == 0 & X4xpff == 0 & X4xpff.but == 0)
ens4 <- subset(ens2,
                 X6xpff.but == 0 & X6xpff.lps == 0 & X6xpff.but.lps == 0)
ens2.nuc <- subset(ens, 
                 quality.removal == 0 & sat.qual.test == 0 & 
                 pff.qual == 0 & but.qual == 0 &
                 attempt != 4 & attempt != 14 & attempt != 6)
ens3.nuc <- subset(ens2.nuc,
                 X6xpff == 0 & X6xpff.but == 0 & X6xpff.lps == 0 & X6xpff.but.lps == 0 & 
                 X2xpff == 0 & X4xpff == 0 & X4xpff.but == 0)
ens4.nuc <- subset(ens2.nuc,
                 X6xpff.but == 0 & X6xpff.lps == 0 & X6xpff.but.lps == 0)


cort2 <- subset(cort0, quality.removal == 0 & sample.size == 0)
cort2.nuc <- subset(cort, quality.removal == 0 & sample.size == 0)
cort2.nuc2 <- subset(cort.m2, quality.removal == 0 & sample.size == 0)

ens_tubb <- subset(ens_tubb0,
                 quality.removal == 0 & sat.qual.test == 0 &
                 attempt.no != 4 & attempt.no != 14 & attempt.no != 6 & attempt.no != 19)

cort_tubb <- subset(cort_tubb0, qual.id == 0 & sat.qual == 0)
```


```{r ens area analysis,  warning = FALSE}
# Define sample size
ens4$m.id <- paste(ens4$attempt,"-",ens4$well)
n.ens <- aggregate(data = ens4,
                          attempt ~ Groups.f,
                          function(attempt) length(unique(attempt)))
m.ens <- aggregate(data = ens4,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
i.ens <- aggregate(data = ens4,
                          number.analyzed ~ Groups.f,
                          function(number.analyzed) length(unique(number.analyzed)))

# sum.avg.int.asyn.ens <- aggregate(data = ens4,
#                           asyn_roi_tot_avg_int ~ Groups.f, FUN = mean)
# sum.avg.int.pff.ens <- aggregate(data = ens4,
#                           pff_roi_tot_avg_int ~ Groups.f, FUN = mean)

# Create table of variables
ens.replicates <- cbind(n.ens[,2],m.ens[,2],i.ens[,2],sum.avg.int.asyn.ens[,2],sum.avg.int.pff.ens[,2])
colnames(ens.replicates) <- c("N","m","no.images","asyn avg int","pff avg int")
rownames(ens.replicates) <-m.ens[,1]
ens.replicates

# Define sample size
cort2$m.id <- paste(cort2$attempt,"-",cort2$well)
n.cort <- aggregate(data = cort2,
                          attempt ~ Groups.f,
                          function(attempt) length(unique(attempt)))
m.cort <- aggregate(data = cort2,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
i.cort <- aggregate(data = cort2,
                          number.analyzed ~ Groups.f,
                          function(number.analyzed) length(unique(number.analyzed)))

# sum.avg.int.asyn.cort <- aggregate(data = cort2,
#                           asyn_roi_tot_avg_int ~ Groups.f, FUN = mean)
# sum.avg.int.pff.cort <- aggregate(data = cort2,
#                           pff_roi_tot_avg_int ~ Groups.f, FUN = mean)

# Create table of variables
cort.replicates <- cbind(n.cort[,2],m.cort[,2],i.cort[,2],sum.avg.int.asyn.cort[,2],sum.avg.int.pff.cort[,2])
colnames(cort.replicates) <- c("N","m","no.images","asyn avg int","pff avg int")
rownames(cort.replicates) <-m.cort[,1]
cort.replicates

# Non-parametric test of pff area by total image
mod.ens.pff.p.area <- kruskal.test(pff_p_area_tot ~ Groups.f, data = ens4)
mod.ens.pff.p.area
ph.ens.pff.p.area <- pairwise.wilcox.test(ens4$pff_p_area_tot, ens4$Groups.f,p.adjust.method = "BH")
ph.ens.pff.p.area 

# Plot values
plot.ens.pff.area <- ggplot(ens4, aes(x=Groups.f, y=pff_p_area_tot, fill=Groups.f)) +
  ylab("Fraction of neural area with elevated PFF tag signal (ATTO-590)") + 
  xlab("Groups") + 
  labs(
    title = "Neuron PFF Area") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 0.5, tip_length = 0.05, vjust=0.4) +
    geom_signif(comparisons=list(c("2XPFF", "Control")), annotations="***",
              y_position = 0.65, tip_length = 0.05, vjust=0.4) +
    geom_signif(comparisons=list(c("4XPFF", "Control")), annotations="***",
              y_position = 0.75, tip_length = 0.05, vjust=0.4) +
    geom_signif(comparisons=list(c("6XPFF", "Control")), annotations="***",
              y_position = 0.95, tip_length = 0.05, vjust=0.4) +
    geom_signif(comparisons=list(c("PFF", "PFF+But")), annotations="0.076",
              y_position = 0.5, tip_length = 0.05, vjust=0) +
    geom_signif(comparisons=list(c("PFF", "PFF+LPS")), annotations="*",
              y_position = 0.58, tip_length = 0.05, vjust=0.4) +
    geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="**",
              y_position = 0.85, tip_length = 0.05, vjust=0.4) +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_dotplot(data = ens4, aes(x = Groups.f, y = pff_p_area_tot),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot.ens.pff.area

mod.ens.asyn.p.area <- kruskal.test(asyn_p_area_tot ~ Groups.f, data = ens4)
ph.ens.asyn.p.area <- pairwise.wilcox.test(ens4$asyn_p_area_tot, ens4$Groups.f,p.adjust.method = "BH")
ph.ens.asyn.p.area

plot.ens.asyn.area <- ggplot(ens4, aes(x=Groups.f, y=asyn_p_area_tot, fill=Groups.f)) +
  ylab("Fraction of neural area with elevated asyn") + 
  xlab("Groups") + 
  labs(
    title = "Neuron a-Syn Area") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 50, tip_length = 0.05, vjust=0.4) +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_dotplot(data = ens4, aes(x = Groups.f, y = asyn_p_area_tot),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.001,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none")

plot.ens.asyn.area + coord_cartesian(ylim = c(0,0.01))
```


```{r ens area corrected for dapi,  warning = FALSE}
# Correct area by nuclei count
ens4.nuc$asyn_p_area_tot.corr <-ens4.nuc$asyn_p_area_tot/ens4.nuc$nuclei.no
ens4.nuc$pff_p_area_tot.corr <-ens4.nuc$pff_p_area_tot/ens4.nuc$nuclei.no

# Non-parametric analysis of area by nuclei count (PFF)
mod.ens.asyn.p.area.corr <- kruskal.test(asyn_p_area_tot.corr ~ Groups.f, data = ens4.nuc)
ph.ens.asyn.p.area.corr <- pairwise.wilcox.test(ens4.nuc$asyn_p_area_tot.corr, ens4.nuc$Groups.f,p.adjust.method = "BH")
ph.ens.asyn.p.area.corr 

# Non-parametric analysis of area by nuclei count (a-Syn)
mod.ens.pff.p.area.corr <- kruskal.test(pff_p_area_tot.corr ~ Groups.f, data = ens4.nuc)
ph.ens.pff.p.area.corr <- pairwise.wilcox.test(ens4.nuc$pff_p_area_tot.corr, ens4.nuc$Groups.f,p.adjust.method = "BH")
ph.ens.pff.p.area.corr 

# Plot values
plot.ens.asyn.area <- ggplot(ens4.nuc, aes(x=Groups.f, y=sqrt(asyn_p_area_tot.corr), fill=Groups.f)) +
  ylab("Sqrt Fraction Area per Nuclei (a-Syn)") + 
  xlab("Groups") + 
  labs(
    title = "a-Syn Area per Nuclei (ENS)") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) +
  scale_fill_viridis_d()

plot.ens.asyn.area

plot.ens.pff.area <- ggplot(ens4.nuc, aes(x=Groups.f, y=sqrt(pff_p_area_tot.corr), fill=Groups.f)) +
  ylab("Sqrt Fraction Area per Nuclei (PFF tag)") + 
  xlab("Groups") + 
  labs(
    title = "PFF Area per Nuclei (ENS)") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 0.05, tip_length = 0.02, size = 0.8, textsize = 8, vjust=0.4) +
    geom_signif(comparisons=list(c("2XPFF", "PFF")), annotations="***",
              y_position = 0.125, tip_length = 0.02, size = 0.8, textsize = 8, vjust=0.4) +
    geom_signif(comparisons=list(c("4XPFF", "2XPFF")), annotations="***",
              y_position = 0.20, tip_length = 0.02, size = 0.8, textsize = 8, vjust=0.4) +
    geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="**",
              y_position = 0.20, tip_length = 0.02, size = 0.8, textsize = 8, vjust=0.4) +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) +
  coord_cartesian(ylim = c(0,0.28)) +
  scale_fill_viridis_d()

plot.ens.pff.area

# Analysis
# To determine whether adjusting for random effects would be valuable, a likelihood ratio test is conducted between models with and without random intercept effects. The simplest test (if any) that is significantly improved from a model without adjustments.

# Statistical tests
mod.ens.asyn.nuc.int <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + nuclei.no + 
                               (1|attempt.no), 
                             data = ens4.nuc)
mod.ens.asyn.nuc.int2 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + nuclei.no + 
                               (1|well), 
                             data = ens4.nuc)
mod.ens.asyn.nuc.int3 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + nuclei.no + 
                               (1|attempt.no) + (1|attempt.no:well), 
                             data = ens4.nuc)
mod.ens.asyn.nuc.int4 <- lm(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + nuclei.no, 
                             data = ens4.nuc)

# Likelihood ratio tests
anova(mod.ens.asyn.nuc.int,mod.ens.asyn.nuc.int2,mod.ens.asyn.nuc.int3,mod.ens.asyn.nuc.int4)
anova(mod.ens.asyn.nuc.int,mod.ens.asyn.nuc.int3)

# Result summary
summary(mod.ens.asyn.nuc.int)

# Run finalized model with transformation prepped for back-transformation
fixed.avg.int.tran <- make.tran("power", (1/6))
warp.asyn.avg.int.nuc <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_roi_tot_avg_int) ~ Groups.f + nuclei.no + 
        (1|attempt.no)
         , data = ens4.nuc))

# Statistical summary
summary(em.ens.asyn.avg.int.nuc <- emmeans(warp.asyn.avg.int.nuc, ~Groups.f, type = "response"))
pwpm(em.ens.asyn.avg.int.nuc)

em.ens.nuc.asyn.tot.int<- plot(em.ens.asyn.avg.int.nuc, type = "response", comparisons = TRUE, plotIT=FALSE) +
  labs(title = "Average Intensity adjusted by nuclei (a-Syn, ENS)",
       x = "Estimated Marginal Mean (average intensity a-Syn)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())
em.ens.nuc.asyn.tot.int

# Statistical tests
mod.ens.pff.nuc.int0 <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no), 
                             data = ens4.nuc)

mod.ens.pff.nuc.int1 <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + nuclei.no + 
                               (1|attempt.no), 
                             data = ens4.nuc)

mod.ens.pff.nuc.int2 <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + nuclei.no + 
                               (1|well), 
                             data = ens4.nuc)

mod.ens.pff.nuc.int3 <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + nuclei.no + 
                               (1|attempt.no) + (1|attempt.no:well), 
                             data = ens4.nuc)

mod.ens.pff.nuc.int4 <- lm(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + nuclei.no, 
                             data = ens4.nuc)

# Likelihood ratio tests
anova(mod.ens.pff.nuc.int0,mod.ens.pff.nuc.int1,mod.ens.pff.nuc.int2,mod.ens.pff.nuc.int3,mod.ens.pff.nuc.int4)
anova(mod.ens.pff.nuc.int1,mod.ens.pff.nuc.int3)
summary(mod.ens.pff.nuc.int1)

# Run finalized model with transformation prepped for back-transformation
fixed.avg.int.tran <- make.tran("power", (1/6))
warp.pff.avg.int.nuc <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_roi_tot_avg_int) ~ Groups.f + nuclei.no + 
        (1|attempt.no)
         , data = ens4.nuc))

# Statistical summary
summary(em.ens.pff.avg.int.nuc <- emmeans(warp.pff.avg.int.nuc, ~Groups.f, type = "response"))
pwpm(em.ens.pff.avg.int.nuc)

# Emmeans plot
em.ens.nuc.pff.tot.int<- plot(em.ens.pff.avg.int.nuc, type = "response", comparisons = TRUE, plotIT=FALSE) +
  labs(title = "Average Intensity adjusted by nuclei (PFF, ENS)",
       x = "Estimated Marginal Mean (average intensity PFF)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())
em.ens.nuc.pff.tot.int
```


```{r ens transformation and analysis avg int,  warning = FALSE}
#Identify whether transformation is needed
#skewness(ens4$asyn_roi_tot_avg_int, na.rm = TRUE)
hist(ens4$asyn_roi_tot_avg_int)
hist(ens4$t.asyn_roi_tot_avg_int <- 
       (ens4$asyn_roi_tot_avg_int)^(1/6))
skewness(ens4$t.asyn_roi_tot_avg_int, na.rm = TRUE)

#Analyze a-Syn
mod.ens.asyn.tot.int <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no), 
                             data = ens4)

mod.ens.asyn.tot.int2 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|well), 
                             data = ens4)

mod.ens.asyn.tot.int3 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no) + (1|attempt.no:well), 
                             data = ens4)

mod.ens.asyn.tot.int4 <- lm(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f, 
                             data = ens4)

# Likelihood ratio tests
anova(mod.ens.asyn.tot.int,mod.ens.asyn.tot.int2,mod.ens.asyn.tot.int3,mod.ens.asyn.tot.int4)
anova(mod.ens.asyn.tot.int,mod.ens.asyn.tot.int3)
summary(mod.ens.asyn.tot.int3)

fixed.avg.int.tran <- make.tran("power", (1/6))
warp.asyn.avg.int.tot <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_roi_tot_avg_int) ~ Groups.f + 
        (1|attempt.no) + (1|attempt.no:well)
         , data = ens4))
summary(em.ens.asyn.avg.int.tot <- emmeans(warp.asyn.avg.int.tot, ~Groups.f, type = "response"))
pwpm(em.ens.asyn.avg.int.tot)
plot(em.ens.asyn.avg.int.tot, comparisons=TRUE, type = "response")
#eff_size(em.ens.asyn.tot.int, sigma = sigma(mod.ens.asyn.tot.int), edf = 11)

#Analyze PFF
mod.ens.pff.tot.int <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no), 
                             data = ens4)

mod.ens.pff.tot.int2 <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|well), 
                             data = ens4)

mod.ens.pff.tot.int3 <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no) + (1|attempt.no:well), 
                             data = ens4)

mod.ens.pff.tot.int4 <- lm(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f, 
                             data = ens4)

# Likelihood ratio tests
anova(mod.ens.pff.tot.int,mod.ens.pff.tot.int2,mod.ens.pff.tot.int3,mod.ens.pff.tot.int4)
anova(mod.ens.pff.tot.int,mod.ens.pff.tot.int3)
summary(mod.ens.pff.tot.int)

fixed.avg.int.tran <- make.tran("power", (1/6))
warp.pff.avg.int.tot <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_roi_tot_avg_int) ~ Groups.f + 
        (1|attempt.no)
         , data = ens4))
summary(em.ens.pff.avg.int.tot <- emmeans(warp.pff.avg.int.tot, ~Groups.f, type = "response"))
pwpm(em.ens.pff.avg.int.tot)
plot(em.ens.pff.avg.int.tot, comparisons=TRUE, type = "response")
#eff_size(em.ens.pff.tot.int, sigma = sigma(mod.ens.pff.tot.int), edf = 11)

# Emmeans plot
em.ens.asyn.tot.int<- plot(em.ens.asyn.avg.int.tot, type = "response", comparisons = TRUE, plotIT=FALSE) +
  labs(title = "Average a-Syn intensity within neuron (ENS)",
       x = "Estimated Marginal Mean (Average Intensity a-Syn)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())
em.ens.asyn.tot.int

em.ens.pff.tot.int<- plot(em.ens.pff.avg.int.tot, type = "response", comparisons = TRUE, plotIT=FALSE) +
  labs(title = "Average PFF intensity within neuron (ENS)",
       x = "Estimated Marginal Mean (Average Intensity PFF)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())
em.ens.pff.tot.int

# Plot values
plot.ens.asyn.int.group <- ggplot(ens4, aes(x=Groups.f, y=asyn_roi_tot_avg_int, fill=attempt.no)) +
  ylab("Average a-Syn Intensity within Neuron") + 
  xlab("Groups") + 
  labs(
    title = "Alpha synuclein average intensity ENS", 
    subtitle = "3-8 replicates, 8-21 wells, 40-90 images, Wilcox Non-parametric Test") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_dotplot(data = ens4, aes(x = Groups.f, y = asyn_roi_tot_avg_int,
                                groups = interaction(Groups.f, attempt.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1)

plot.ens.asyn.int.group

ens3.simp <- subset(ens3, (Groups.f == "Control" | Groups.f == "PFF" | Groups.f == "But" | Groups.f == "PFF+But") )

plot.ens.pff.int.group <- ggplot(ens3.simp, aes(x=Groups.f, y=pff_roi_tot_avg_int, fill=attempt.no)) +
  ylab("Average PFF Intensity within Neuron") + 
  xlab("Groups") + 
  labs(
    title = "Variation exists within group and experimental replicate") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_dotplot(data = ens3.simp, aes(x = Groups.f, y = pff_roi_tot_avg_int,
                                groups = interaction(Groups.f, attempt.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1)

plot.ens.pff.int.group

plot.ens.asyn.int.group.simp <- ggplot(ens3, aes(x=Groups.f, y=asyn_roi_tot_avg_int, fill=attempt.no)) +
  ylab("Average a-Syn Intensity within Neuron") + 
  xlab("Groups") + 
  labs(
    title = "Neuron a-Syn Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_dotplot(data = ens3, aes(x = Groups.f, y = asyn_roi_tot_avg_int,
                                groups = interaction(Groups.f, attempt.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1)

plot.ens.asyn.int.group.simp
```


```{r nuclei count by group ens,  warning = FALSE}
# Plot values
plot.ens.nuclei <- ggplot(ens4.nuc, aes(x=Groups.f, y=nuclei.no, fill=Groups.f)) +
  ylab("Number of Nuclei per image") + 
  xlab("Groups") + 
  labs(
    title = "Comparison of Nuclei Count per Groups")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_dotplot(data = ens4.nuc, aes(x = Groups.f, y = nuclei.no),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none", text = element_text(size = 14))

plot.ens.nuclei


plot.ens.nuclei2 <- ggplot(ens4.nuc, aes(x=Groups.f, y=asyn_roi_tot_avg_int, fill=bin.nuclei)) +
  ylab("Average Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Comparison of Nuclei Count per Groups")+
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_dotplot(data = ens4.nuc, aes(x = Groups.f, y = asyn_roi_tot_avg_int,
                                groups = interaction(Groups.f, nuclei.no)),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(text = element_text(size = 14))

plot.ens.nuclei2

plot.ens.nuclei3 <- ggplot(ens4.nuc, aes(x=nuclei.no, y=asyn_roi_tot_avg_int, colour=Groups.f)) +
  geom_point() +
  xlab("Number of Nuclei in image") +
  ylab("Average a-Syn Intensity in Neuron") +
  theme(text = element_text(size = 14))

plot.ens.nuclei3

ens4.nuc.2 <- subset(ens4.nuc, 
                 (Groups.f == "Control" | Groups.f == "But" | Groups.f == "PFF" | Groups.f == "PFF+But"))

ens4.nuc.3 <- subset(ens4.nuc, 
                 Groups.f == "Control")

plot.ens.nuclei3.2 <- ggplot(ens4.nuc.2, aes(x=nuclei.no, y=asyn_roi_tot_avg_int, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Nuclei count vs. a-Syn intensity within Neuron (all groups)"
  ) +
  xlab("Number of Nuclei in image") +
  ylab("Average a-Syn Intensity in Neuron") +
  theme(text = element_text(size = 14))

plot.ens.nuclei3.2

plot.ens.nuclei3.2p <- ggplot(ens4.nuc.2, aes(x=nuclei.no, y=pff_roi_tot_avg_int, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Nuclei count vs. PFF intensity within Neuron (all groups)",
    col = "Groups") +
  xlab("Number of Nuclei in image") +
  ylab("Average PFF Intensity in Neuron") +
  theme(text = element_text(size = 14))

plot.ens.nuclei3.2p

plot.ens.comp3.3 <- ggplot(ens4.nuc.2, aes(x=pff_roi_tot_avg_int, y=asyn_roi_tot_avg_int, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Comparison of a-Syn intensity with PFF intensity",
    col = "Groups") +
  xlab("Average PFF Intensity in Neuron") +
  ylab("Average a-Syn Intensity in Neuron") +
  theme(text = element_text(size = 14))

plot.ens.comp3.3

#Scatterplot Sample Size
ens4.nuc.2$m.id <- paste(ens4.nuc.2$attempt,"-",ens4.nuc.2$well)
n.ens.nuc <- aggregate(data = ens4.nuc.2,
                          attempt ~ Groups.f,
                          function(attempt) length(unique(attempt)))
m.ens.nuc <- aggregate(data = ens4.nuc.2,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
i.ens.nuc <- aggregate(data = ens4.nuc.2,
                          number.analyzed ~ Groups.f,
                          function(number.analyzed) length(unique(number.analyzed)))
ens.nuc.replicates <- cbind(n.ens.nuc[,2],m.ens.nuc[,2],i.ens.nuc[,2])
colnames(ens.nuc.replicates) <- c("N","m","no.images")
rownames(ens.nuc.replicates) <- c("Control","But","PFF","PFF+But")
ens.nuc.replicates

#Corrected Plots Sample Size
ens4.nuc$m.id <- paste(ens4.nuc$attempt,"-",ens4.nuc$well)
n.ens.nuc <- aggregate(data = ens4.nuc,
                          attempt ~ Groups.f,
                          function(attempt) length(unique(attempt)))
m.ens.nuc <- aggregate(data = ens4.nuc,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
i.ens.nuc <- aggregate(data = ens4.nuc,
                          number.analyzed ~ Groups.f,
                          function(number.analyzed) length(unique(number.analyzed)))
ens.nuc.replicates.corr <- cbind(n.ens.nuc[,2],m.ens.nuc[,2],i.ens.nuc[,2])
colnames(ens.nuc.replicates.corr) <- c("N","m","no.images")
rownames(ens.nuc.replicates.corr) <- n.ens.nuc[,1]
ens.nuc.replicates.corr

plot.ens.nuclei3.3p <- ggplot(ens4.nuc.3, aes(x=nuclei.no, y=pff_roi_tot_avg_int, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Nuclei count vs. PFF intensity within Neuron",
    col = "Groups") +
  xlab("Number of Nuclei in image") +
  ylab("Average PFF Intensity in Neuron") +
  theme(legend.position="none", text = element_text(size = 14))

plot.ens.nuclei3.3p

plot.ens.nuclei3.3a <- ggplot(ens4.nuc.3, aes(x=nuclei.no, y=asyn_roi_tot_avg_int, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Nuclei count vs. a-Syn intensity within Neuron",
    col = "Groups") +
  xlab("Number of Nuclei in image") +
  ylab("Average a-Syn Intensity within Neuron") +
  theme(legend.position="none", text = element_text(size = 14))

plot.ens.nuclei3.3a

plot.ens.nuclei3.4p <- ggplot(ens4.nuc.3, aes(x=nuclei.no, y=pff_p_area_tot, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Nuclei count vs. PFF elevated area",
    col = "Groups") +
  xlab("Number of Nuclei in image") +
  ylab("Fraction of neural area with elevated PFF (ATTO-590)") +
  theme(legend.position="none",text = element_text(size = 14))

plot.ens.nuclei3.4p

plot.ens.nuclei3.4a <- ggplot(ens4.nuc.3, aes(x=nuclei.no, y=asyn_p_area_tot, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Nuclei count vs. a-Syn elevated area",
    col = "Groups") +
  xlab("Number of Nuclei in image") +
  ylab("Fraction of neural area with elevated a-Syn") +
  theme(legend.position="none",text = element_text(size = 14))

plot.ens.nuclei3.4a

plot.ens.nuclei3.5p <- ggplot(ens4.nuc.2, aes(x=nuclei.no, y=pff_p_area_tot, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Nuclei count vs. fraction elevated PFF area with Neuron (all groups)",
    col = "Groups") +
  xlab("Number of Nuclei in image") +
  ylab("Average PFF elevated area in Neuron") +
  theme(text = element_text(size = 14))

plot.ens.nuclei3.5p

plot.ens.nuclei3.5a <- ggplot(ens4.nuc.2, aes(x=nuclei.no, y=asyn_p_area_tot, colour=Groups.f)) +
  geom_point() +
  labs(
    title = "Nuclei count vs. fraction elevated a-Syn area with Neuron (all groups)",
    col = "Groups") +
  xlab("Number of Nuclei in image") +
  ylab("Average a-Syn elevated area in Neuron") +
  theme(text = element_text(size = 14))

plot.ens.nuclei3.5a


mod.ens.nuclei.no3a <- lm(asyn_roi_tot_avg_int^(1/6) ~ nuclei.no,data = ens4.nuc.3)
summary(mod.ens.nuclei.no3a)


mod.ens.nuclei.no3p <- lm(pff_roi_tot_avg_int^(1/6) ~ nuclei.no,data = ens4.nuc.3)
summary(mod.ens.nuclei.no3p)

mod.ens.nuclei.no2b <- lm(asyn_roi_tot_avg_int^(1/6) ~ bin.nuclei * Groups.f,data = ens4.nuc.2)
em.ens.nuclei.no2b <- emmeans(mod.ens.nuclei.no2b, specs = "Groups.f", by="bin.nuclei")
pwpm(em.ens.nuclei.no2b)


mod.ens.nuc.asyn.tot.int <- lm(asyn_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no,data = ens4.nuc)
summary(mod.ens.nuclei.no3a)


mod.ens.nuclei.no3p <- lm(pff_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no,data = ens4.nuc)
summary(mod.ens.nuclei.no3p)
```

```{r nuclei adjusted,  warning = FALSE}
#Subset
ens5.nuc <- subset(ens4.nuc, nuclei.no < 20)

#Avg Int
mod.asyn.ens.nuc.int <- lmer(asyn_roi_line_avg_int^(1/6) ~ 1 + Groups.f + (1|attempt.no), data = ens5.nuc)
warp.ens.asyn.avg.int.raw <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_roi_line_avg_int) ~ Groups.f + 
        (1|attempt.no)
         , data = ens5.nuc))
summary(em.asyn.ens.nuc.int <- emmeans(warp.ens.asyn.avg.int.raw, ~Groups.f, type = "response"))
pwpm(em.asyn.ens.nuc.int)

mod.pff.ens.nuc.int <- lmer(pff_roi_line_avg_int^(1/6) ~ 1 + Groups.f + (1|attempt.no), data = ens5.nuc)
warp.ens.pff.avg.int.raw <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_roi_line_avg_int) ~ Groups.f + 
        (1|attempt.no)
         , data = ens5.nuc))
summary(em.pff.ens.nuc.int <- emmeans(warp.ens.pff.avg.int.raw, ~Groups.f, type = "response"))
pwpm(em.pff.ens.nuc.int)

#Area
mod.ens5.nuc.asyn.p.area.corr <- kruskal.test(asyn_p_area_tot.corr ~ Groups.f, data = ens5.nuc)
ph.ens5.nuc.asyn.p.area.corr <- pairwise.wilcox.test(ens5.nuc$asyn_p_area_tot.corr, ens5.nuc$Groups.f,p.adjust.method = "BH")
ph.ens5.nuc.asyn.p.area.corr 

mod.ens5.nuc.pff.p.area <- kruskal.test(pff_p_area_tot ~ Groups.f, data = ens5.nuc)
ph.ens5.nuc.pff.p.area <- pairwise.wilcox.test(ens5.nuc$pff_p_area_tot, ens5.nuc$Groups.f,p.adjust.method = "BH")
ph.ens5.nuc.pff.p.area 

```


```{r cortical,  warning = FALSE}
#Sample Size
cort2$m.id <- paste(cort2$attempt,"-",cort2$well)
n.cort <- aggregate(data = cort2,
                          attempt ~ Groups.f,
                          function(attempt) length(unique(attempt)))
m.cort <- aggregate(data = cort2,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
i.cort <- aggregate(data = cort2,
                          number.analyzed ~ Groups.f,
                          function(number.analyzed) length(unique(number.analyzed)))
cort.replicates <- cbind(n.cort[,2],m.cort[,2],i.cort[,2])
colnames(cort.replicates) <- c("N","m","no.images")
rownames(cort.replicates) <- c("Control","PFF","4XPFF")
cort.replicates

#Transformation
hist(cort2$pff_roi_tot_avg_int)
hist(cort2$t.pff_roi_tot_avg_int <-
       (cort2$pff_roi_tot_avg_int)^(1/6))

#Analyze a-Syn
mod.cor.asyn.tot.int <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no), 
                             data = cort2)

mod.cor.asyn.tot.int2 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|well), 
                             data = cort2)

mod.cor.asyn.tot.int3 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no) + (1|attempt.no:well), 
                             data = cort2)

mod.cor.asyn.tot.int4 <- lm(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f, 
                             data = cort2)

anova(mod.cor.asyn.tot.int,mod.cor.asyn.tot.int2,mod.cor.asyn.tot.int3,mod.cor.asyn.tot.int4)
anova(mod.cor.asyn.tot.int,mod.cor.asyn.tot.int3)
summary(mod.cor.asyn.tot.int)

fixed.avg.int.tran <- make.tran("power", (1/6))
warp.cort.asyn.avg.int.tot <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_roi_tot_avg_int) ~ Groups.f + 
        (1|attempt.no)
         , data = cort2))
summary(em.cor.asyn.avg.int.tot <- emmeans(warp.cort.asyn.avg.int.tot, ~Groups.f, type = "response"))
pwpm(em.cor.asyn.avg.int.tot)
plot(em.cor.asyn.avg.int.tot, comparisons=TRUE, type = "response")
#eff_size(em.cor.asyn.avg.int.tot, sigma = sigma(mod.cor.asyn.tot.int), edf = 2)

#Analyze PFF
mod.cor.pff.tot.int <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no), 
                             data = cort2)

mod.cor.pff.tot.int2 <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|well), 
                             data = cort2)

mod.cor.pff.tot.int3 <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                               (1|attempt.no) + (1|attempt.no:well), 
                             data = cort2)

mod.cor.pff.tot.int4 <- lm(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f, 
                             data = cort2)

anova(mod.cor.pff.tot.int,mod.cor.pff.tot.int2,mod.cor.pff.tot.int3,mod.cor.pff.tot.int4)
anova(mod.cor.pff.tot.int,mod.cor.pff.tot.int3)
summary(mod.cor.pff.tot.int)

fixed.avg.int.tran <- make.tran("power", (1/6))
warp.pff.avg.int.tot <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_roi_tot_avg_int) ~ Groups.f + 
        (1|attempt.no)
         , data = cort2))
summary(em.cor.pff.avg.int.tot <- emmeans(warp.pff.avg.int.tot, ~Groups.f, type = "response"))
pwpm(em.cor.pff.avg.int.tot)
plot(em.cor.pff.avg.int.tot, comparisons=TRUE, type = "response")
#eff_size(em.cor.pff.tot.int, sigma = sigma(mod.cor.pff.tot.int), edf = 2)

em.cor.asyn.tot.int<- plot(em.cor.asyn.avg.int.tot, type = "response", comparisons = TRUE, plotIT=FALSE) +
  labs(title = "Average a-Syn intensity within neuron (Cortex)",
       x = "Estimated Marginal Mean (average intensity a-Syn)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())
em.cor.asyn.tot.int

em.cor.pff.tot.int<- plot(em.cor.pff.avg.int.tot, type = "response", comparisons = TRUE, plotIT=FALSE) +
  labs(title = "Average PFF tag intensity within neuron (Cortex)",
       x = "Estimated Marginal Mean (average intensity PFF)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())
em.cor.pff.tot.int

plot.cort.asyn.int <- ggplot(cort2, aes(x=Groups.f, y=asyn_roi_tot_avg_int)) +
  ylab("Average a-Syn Intensity within Neuron") + 
  xlab("Groups") + 
  labs(
    title = "Average a-Syn intensity within neuron (Cortex)") +
  scale_x_discrete(labels=c("Control", "1 µg PFF", "4µg PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean',fill = "aquamarine4") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 75, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "Control")), annotations="***",
              y_position = 82, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "4XPFF")), annotations="**",
              y_position = 75, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = cort2, aes(x = Groups.f, y = asyn_roi_tot_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.6,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.title.x = element_blank())

plot.cort.asyn.int

plot.cort.pff.int <- ggplot(cort2, aes(x=Groups.f, y=pff_roi_tot_avg_int)) +
  ylab("Average PFF (ATTO-590) Intensity within Neuron") + 
  xlab("Groups") + 
  labs(
    title = "Cortical Average Intensity with Neuron (PFF)") +
  scale_x_discrete(labels=c("Control", "1 µg PFF", "4µg PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean',fill = "aquamarine4") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 80, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "Control")), annotations="***",
              y_position = 87, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "4XPFF")), annotations="***",
              y_position = 80, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = cort2, aes(x = Groups.f, y = pff_roi_tot_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.35,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.title.x = element_blank())

plot.cort.pff.int
```

```{r cort area,  warning = FALSE}
# Plot histogram
hist(cort2$asyn_p_area_tot)
hist(cort2$t.asyn_p_area_tot <-
       1/(cort2$asyn_p_area_tot+0.00001))
skewness(cort2$t.asyn_p_area_tot, na.rm = FALSE)

# Non-parametric test
mod.cor.asyn.p.area <- kruskal.test(asyn_p_area_tot ~ Groups.f, data = cort2)
ph.cor.asyn.p.area <- pairwise.wilcox.test(cort2$asyn_p_area_tot, cort2$Groups.f,p.adjust.method = "BH")
ph.cor.asyn.p.area 

# Non-parametric test
mod.cor.pff.p.area <- kruskal.test(pff_p_area_tot ~ Groups.f, data = cort2)
ph.cor.pff.p.area <- pairwise.wilcox.test(cort2$pff_p_area_tot, cort2$Groups.f,p.adjust.method = "BH")
ph.cor.pff.p.area 

# Plot values
plot.cort.asyn.area <- ggplot(cort2, aes(x=Groups.f, y=sqrt(asyn_p_area_tot))) +
  ylab("Sqrt Fraction Area within Neuron (a-Syn)") + 
  xlab("Groups") + 
  labs(
    title = "Cortex a-Syn Area")+
  scale_x_discrete(labels=c("Control", "1 µg PFF", "4µg PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean', fill = "aquamarine4") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black")  +
  geom_signif(comparisons=list(c("4XPFF", "Control")), annotations="***",
              y_position = 0.2, tip_length = 0.01, size = 0.8, textsize = 8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "4XPFF")), annotations="***",
              y_position = 0.23, tip_length = 0.01, size = 0.8, textsize = 8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.title.x = element_blank())

plot.cort.asyn.area 

# cort3 <- subset(cort2, attempt.no != 17)
# plot.cort.asyn.area <- ggplot(cort3, aes(x=Groups.f, y=asyn_p_area_tot, fill=Groups.f)) +
#   scale_x_discrete(labels=c("Control", "1 µg PFF", "4µg PFF")) +
#   geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
#     stat_summary(fun.data = mean_se, geom = "errorbar", position = 'dodge', colour = "black")
# 
# plot.cort.asyn.area

plot.cort.pff.area <- ggplot(cort2, aes(x=Groups.f, y=sqrt(pff_p_area_tot))) +
  ylab("Sqrt Fraction Area within Neuron (PFF)") + 
  xlab("Groups") + 
  labs(
    title = "Cortex PFF Area")+
  scale_x_discrete(labels=c("Control", "1 µg PFF", "4µg PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean', fill = "aquamarine4") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 0.625, tip_length = 0.05, size = 0.8, textsize = 8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "Control")), annotations="***",
              y_position = 0.725, tip_length = 0.05, size = 0.8, textsize = 8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "4XPFF")), annotations="***",
              y_position = 0.625, tip_length = 0.05, size = 0.8, textsize = 8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.title.x = element_blank()) 

plot.cort.pff.area

```

```{r cort avg int corrected for nuclei}
#a-Syn models
mod.cor.asyn.tot.int.corr <- lm(asyn_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no, data = cort2.nuc)
mod.cor.asyn.tot.int.corr2 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no +
                                     (1|attempt.no),
                                   data = cort2.nuc)
mod.cor.asyn.tot.int.corr3 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no +
                                     (1|well),
                                   data = cort2.nuc)
mod.cor.asyn.tot.int.corr4 <- lmer(asyn_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no +
                                     (1|attempt.no) + (1|attempt:well),
                                   data = cort2.nuc)

# Assess best models
lrtest(mod.cor.asyn.tot.int.corr,mod.cor.asyn.tot.int.corr2,
      mod.cor.asyn.tot.int.corr3,mod.cor.asyn.tot.int.corr4)
anova(mod.cor.asyn.tot.int.corr2,mod.cor.asyn.tot.int.corr3,mod.cor.asyn.tot.int.corr4)

# Run optimal models set for back-transformation
warp.cor.asyn.tot.int.corr <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_roi_tot_avg_int) ~ Groups.f + nuclei.no +
        (1|attempt.no) + (1|attempt:well)
         , data = cort2.nuc))

# Statistical Results
summary(cor.asyn.tot.int.corr <- emmeans(warp.cor.asyn.tot.int.corr, ~Groups.f, type = "response"))
pwpm(cor.asyn.tot.int.corr)

# Emmeans
em.cor.asyn.tot.int.corr <- emmeans(mod.cor.asyn.tot.int.corr, ~"Groups.f", )
pwpm(em.cor.asyn.tot.int.corr)

# Plot emmeans
em.cor.asyn.tot.int.corr<- plot(em.cor.asyn.tot.int.corr, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Average Intensity adjusted by nuclei (a-Syn, Cortex)",
       x = "Estimated Marginal Mean (Average Intensity a-Syn)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())
em.cor.asyn.tot.int.corr

#PFF models
mod.cor.pff.tot.int.corr <- lm(pff_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no, data = cort2.nuc)
mod.cor.pff.tot.int.corr2 <- lmer(pff_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no +
                                     (1|attempt.no),
                                   data = cort2.nuc)
mod.cor.pff.tot.int.corr3 <- lmer(pff_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no +
                                     (1|well),
                                   data = cort2.nuc)
mod.cor.pff.tot.int.corr4 <- lmer(pff_roi_tot_avg_int^(1/6) ~ Groups.f + nuclei.no +
                                     (1|attempt.no) + (1|attempt:well),
                                   data = cort2.nuc)

# Assess best models
lrtest(mod.cor.pff.tot.int.corr,mod.cor.pff.tot.int.corr2,
      mod.cor.pff.tot.int.corr3,mod.cor.pff.tot.int.corr4)
anova(mod.cor.pff.tot.int.corr2,mod.cor.pff.tot.int.corr3,mod.cor.pff.tot.int.corr4)

# Run optimal models set for back-transformation
warp.cor.pff.tot.int.corr <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_roi_tot_avg_int) ~ Groups.f + nuclei.no +
        (1|attempt.no) + (1|attempt:well)
         , data = cort2.nuc))

# Statistical Results
summary(cor.pff.tot.int.corr <- emmeans(warp.cor.pff.tot.int.corr, ~Groups.f, type = "response"))
pwpm(cor.pff.tot.int.corr)

# Emmeans
em.cor.pff.tot.int.corr <- emmeans(mod.cor.pff.tot.int.corr, ~"Groups.f", )
pwpm(em.cor.pff.tot.int.corr)

# Plot emmeans
em.cor.pff.tot.int.corr<- plot(em.cor.pff.tot.int.corr, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Average Intensity adjusted by nuclei (PFF, Cortex)",
       x = "Estimated Marginal Mean (Average Intensity PFF)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())
em.cor.pff.tot.int.corr

```

```{r cort area corrected for nuclei,  warning = FALSE}
# Calculate nuclei normalized area
cort2.nuc$asyn_p_area_tot.corr <-cort2.nuc$asyn_p_area_tot/cort2.nuc$nuclei.no
cort2.nuc$pff_p_area_tot.corr <-cort2.nuc$pff_p_area_tot/cort2.nuc$nuclei.no

# Non-parametric tests
mod.cor.asyn.p.area.corr <- kruskal.test(asyn_p_area_tot.corr ~ Groups.f, data = cort2.nuc)
ph.cor.asyn.p.area.corr <- pairwise.wilcox.test(
    cort2.nuc$asyn_p_area_tot.corr, 
    cort2.nuc$Groups.f,
    p.adjust.method = "BH")
ph.cor.asyn.p.area.corr 

mod.cor.pff.p.area.corr <- kruskal.test(pff_p_area_tot.corr ~ Groups.f, data = cort2.nuc)
ph.cor.pff.p.area.corr <- pairwise.wilcox.test(
    cort2.nuc$pff_p_area_tot.corr, 
    cort2.nuc$Groups.f,
    p.adjust.method = "BH")
ph.cor.pff.p.area.corr 

# Plot values
plot.cort.asyn.area <- ggplot(cort2.nuc, aes(x=Groups.f, y=sqrt(asyn_p_area_tot.corr))) +
  ylab("Sqrt Fraction Area with Neuron per nuclei (a-Syn)") + 
  xlab("Groups") + 
  labs(
    title = "a-Syn Area per Nuclei (Cortex)")+
  scale_x_discrete(labels=c("Control", "1 µg PFF", "4µg PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean',fill = "aquamarine4") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="0.06",
              y_position = 0.04, tip_length = 0.005, size = 0.8, textsize = 4, vjust=0) +
  geom_signif(comparisons=list(c("4XPFF", "Control")), annotations="***",
              y_position = 0.05, tip_length = 0.005, size = 0.8, textsize = 8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "4XPFF")), annotations="***",
              y_position = 0.04, tip_length = 0.005, size = 0.8, textsize = 8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,0.09))

plot.cort.asyn.area

# cort3 <- subset(cort2.nuc, attempt.no != 17)
# plot.cort.asyn.area <- ggplot(cort3, aes(x=Groups.f, y=asyn_p_area_tot.corr, fill=Groups.f)) +
#   scale_x_discrete(labels=c("Control", "1 µg PFF", "4µg PFF")) +
#   geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
#     stat_summary(fun.data = mean_se, geom = "errorbar", position = 'dodge', colour = "black")
# 
# plot.cort.asyn.area

plot.cort.pff.area <- ggplot(cort2.nuc, aes(x=Groups.f, y=sqrt(pff_p_area_tot.corr))) +
  ylab("Sqrt Fraction Area with Neuron per nuclei (PFF tag)") + 
  xlab("Groups") + 
  labs(
    title = "PFF Area per Nuclei (Cortex)")+
  scale_x_discrete(labels=c("Control", "1 µg PFF", "4µg PFF")) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean',fill = "aquamarine4") +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 0.30, tip_length = 0.02, size = 0.8, textsize = 8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "Control")), annotations="***",
              y_position = 0.35, tip_length = 0.02, size = 0.8, textsize = 8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "4XPFF")), annotations="***",
              y_position = 0.30, tip_length = 0.02, size = 0.8, textsize = 8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.title.x = element_blank()) 
  
plot.cort.pff.area

# Sample Size Calculation
cort2.nuc$m.id <- paste(cort2.nuc$attempt,"-",cort2.nuc$well)
n.cort.nuc <- aggregate(data = cort2.nuc,
                          attempt ~ Groups.f,
                          function(attempt) length(unique(attempt)))
m.cort.nuc <- aggregate(data = cort2.nuc,
                          m.id ~ Groups.f,
                          function(m.id) length(unique(m.id)))
i.cort.nuc <- aggregate(data = cort2.nuc,
                          number.analyzed ~ Groups.f,
                          function(number.analyzed) length(unique(number.analyzed)))

# Merge Values to a Table
cort.nuc.replicates.corr <- cbind(n.cort.nuc[,2],m.cort.nuc[,2],i.cort.nuc[,2])
colnames(cort.nuc.replicates.corr) <- c("N","m","no.images")
rownames(cort.nuc.replicates.corr) <- n.cort.nuc[,1]
cort.nuc.replicates.corr
```


```{r ens optimization raw}
#Elevated Area ENS non-parateric tests
mod.asyn.p.raw.area <- kruskal.test(asyn_p_area ~ Groups.f, data = ens4)
ph.asyn.p.raw.area <- pairwise.wilcox.test(ens4$asyn_p_area, ens4$Groups.f,p.adjust.method = "BH")
ph.asyn.p.raw.area

mod.pff.p.raw.area  <- kruskal.test(pff_p_area ~ Groups.f, data = ens4)
ph.pff.p.raw.area <- pairwise.wilcox.test(ens4$pff_p_area, ens4$Groups.f,p.adjust.method = "BH")
ph.pff.p.raw.area

# Plot values
plot.ens.asyn.raw.area <- ggplot(ens4, aes(x=Groups.f, y=sqrt(asyn_p_area), fill=Groups.f)) +
  ylab("Sqrt fraction area within whole image (a-Syn)") + 
  xlab("Groups") + 
  labs(
    title = "Whole Image a-Syn Area") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("Control", "4XPFF")), annotations="*",
              y_position = 0.0035, tip_length = 0.002, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "6XPFF")), annotations="*",
              y_position = 0.006, tip_length = 0.002, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="*",
              y_position = 0.001, tip_length = 0.002, size = 0.8, textsize=8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,0.03)) +
  scale_fill_viridis_d()

plot.ens.asyn.raw.area

plot.ens.pff.raw.area <- ggplot(ens4, aes(x=Groups.f, y=sqrt(pff_p_area), fill=Groups.f)) +
  ylab("Sqrt fraction area within whole image (PFF tag)") + 
  xlab("Groups") + 
  labs(
    title = "Whole Image PFF Area") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 0.1, tip_length = 0.01, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "PFF")), annotations="***",
              y_position = 0.15, tip_length = 0.01, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="***",
              y_position = 0.24, tip_length = 0.01, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="***",
              y_position = 0.1, tip_length = 0.01, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="***",
              y_position = 0.24, tip_length = 0.01, size = 0.8, textsize=8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  scale_fill_viridis_d()

plot.ens.pff.raw.area

#Avg Int Raw Models
mod.asyn.raw.int0 <- lm(asyn_avg_int^(1/6) ~ 1 + Groups.f, data = ens4)
mod.asyn.raw.int1 <- lmer(asyn_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no), data = ens4)
mod.asyn.raw.int2 <- lmer(asyn_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|well), data = ens4)
mod.asyn.raw.int3 <- lmer(asyn_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no) + (1|attempt.no:well), data = ens4)

# Identify Best Models
lrtest(mod.asyn.raw.int0,mod.asyn.raw.int1,mod.asyn.raw.int2,mod.asyn.raw.int3)
anova(mod.asyn.raw.int1,mod.asyn.raw.int2,mod.asyn.raw.int3)

mod.asyn.raw.int <- lmer(asyn_avg_int^(1/6) ~ 1 + Groups.f + (1|attempt.no), data = ens4)

# Develop Avg Int a-Syn model designed for back-transformation
warp.ens.asyn.avg.int.raw <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_avg_int) ~ Groups.f + 
        (1|attempt.no) + (1|attempt.no:well)
         , data = ens4))

# Statistical Summary
summary(em.asyn.raw.int <- emmeans(warp.ens.asyn.avg.int.raw, ~Groups.f, type = "response"))
pwpm(em.asyn.raw.int)

# PFF Models
mod.pff.raw.int0 <- lm(pff_avg_int^(1/6) ~ 1 + Groups.f, data = ens4)
mod.pff.raw.int1 <- lmer(pff_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no), data = ens4)
mod.pff.raw.int2 <- lmer(pff_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|well), data = ens4)
mod.pff.raw.int3 <- lmer(pff_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no) + (1|attempt.no:well), data = ens4)

# Identify Optimal Model
lrtest(mod.pff.raw.int0,mod.pff.raw.int1,mod.pff.raw.int2,mod.pff.raw.int3)
anova(mod.pff.raw.int1,mod.pff.raw.int2,mod.pff.raw.int3)

mod.pff.raw.int <- lmer(pff_avg_int^(1/6) ~ 1 + Groups.f + (1|attempt.no), data = ens4)

# Develop Avg Int PFF model designed for back-transformation
warp.ens.pff.avg.int.raw <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_avg_int) ~ Groups.f + 
        (1|attempt.no) + (1|attempt.no:well)
         , data = ens4))

# Statistical Summary
summary(em.pff.raw.int <- emmeans(warp.ens.pff.avg.int.raw, ~Groups.f, type = "response"))
pwpm(em.pff.raw.int)

# Plot Values
plot.ens.asyn.raw.int <- ggplot(ens4, aes(x=Groups.f, y=asyn_avg_int, fill=Groups.f)) +
  ylab("Average a-Syn Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Whole Image a-Syn Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("LPS", "Control")), annotations="***",
              y_position = 28, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("LPS", "But+LPS")), annotations="**",
              y_position = 28, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "PFF")), annotations="***",
              y_position = 32, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "2XPFF")), annotations="***",
              y_position = 36, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "6XPFF")), annotations="***",
              y_position = 40, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = ens4, aes(x = Groups.f, y = asyn_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  coord_cartesian(ylim = c(0,43)) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  scale_fill_viridis_d()

plot.ens.asyn.raw.int

plot.ens.pff.raw.int <- ggplot(ens4, aes(x=Groups.f, y=pff_avg_int, fill=Groups.f)) +
  ylab("Average PFF (ATTO-590) Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Whole Image PFF Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 20, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="***",
              y_position = 20, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="***",
              y_position = 35, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "Control")), annotations="***",
              y_position = 28, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="***",
              y_position = 35, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "6XPFF")), annotations="**",
              y_position = 39, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = ens4, aes(x = Groups.f, y = pff_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,42)) +
  scale_fill_viridis_d()

plot.ens.pff.raw.int

em.ens.asyn.raw.int<- plot(em.asyn.raw.int, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Whole Image a-Syn Average Intensity",
       x = "Estimated Marginal Mean (average intensity a-Syn)",
       y = "Groups")
em.ens.asyn.raw.int

em.ens.pff.raw.int<- plot(em.pff.raw.int, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Whole Image PFF Average Intensity",
       x = "Estimated Marginal Mean (average intensity PFF)",
       y = "Groups")
em.ens.pff.raw.int
```

```{r ens optimization neurite}
# ENS Area Non-parametric tests
mod.asyn.p.line.area <- kruskal.test(asyn_p_area_line ~ Groups.f, data = ens4)
ph.asyn.p.line.area <- pairwise.wilcox.test(ens4$asyn_p_area_line, ens4$Groups.f,p.adjust.method = "BH")
ph.asyn.p.line.area

mod.pff.p.line.area  <- kruskal.test(pff_p_area_line ~ Groups.f, data = ens4)
ph.pff.p.line.area <- pairwise.wilcox.test(ens4$pff_p_area_line, ens4$Groups.f,p.adjust.method = "BH")
ph.pff.p.line.area

# Plot Values
plot.ens.asyn.line.area <- ggplot(ens4, aes(x=Groups.f, y=sqrt(asyn_p_area_line), fill=Groups.f)) +
  ylab("Sqrt fraction area with neurite (a-Syn)") + 
  xlab("Groups") + 
  labs(
    title = "Neurite a-Syn Area") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
  geom_signif(comparisons=list(c("Control", "But")), annotations="*",
              y_position = 0.001, tip_length = 0.002, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "4XPFF")), annotations="*",
              y_position = 0.01, tip_length = 0.002, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "6XPFF")), annotations="*",
              y_position = 0.013, tip_length = 0.002, size = 0.8, textsize=8, vjust=0.4) +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,0.038)) +
  scale_fill_viridis_d()

plot.ens.asyn.line.area

plot.ens.pff.line.area <- ggplot(ens4, aes(x=Groups.f, y=sqrt(pff_p_area_line), fill=Groups.f)) +
  ylab("Sqrt fraction area with neurite (PFF tag)") + 
  xlab("Groups") + 
  labs(
    title = "Neurite PFF Area") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 0.2, tip_length = 0.015, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "PFF")), annotations="*",
              y_position = 0.3, tip_length = 0.015, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="***",
              y_position = 0.48, tip_length = 0.015, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="**",
              y_position = 0.48, tip_length = 0.015, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But")), annotations="0.0546",
              y_position = 0.15, tip_length = 0.015, size = 0.8, textsize=4, vjust=0) +
  geom_signif(comparisons=list(c("PFF", "PFF+LPS")), annotations="*",
              y_position = 0.2, tip_length = 0.015, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="**",
              y_position = 0.25, tip_length = 0.015, size = 0.8, textsize=8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,0.55)) +
  scale_fill_viridis_d()
  
plot.ens.pff.line.area

#Avg Int Neurite Models
mod.asyn.line.int0 <- lm(asyn_roi_line_avg_int^(1/6) ~ 1 + Groups.f, data = ens4)
mod.asyn.line.int1 <- lmer(asyn_roi_line_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no), data = ens4)
mod.asyn.line.int2 <- lmer(asyn_roi_line_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|well), data = ens4)
mod.asyn.line.int3 <- lmer(asyn_roi_line_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no) + (1|attempt.no:well), data = ens4)

# Identify Best Models
lrtest(mod.asyn.line.int0,mod.asyn.line.int1,mod.asyn.line.int2,mod.asyn.line.int3)
anova(mod.asyn.line.int1,mod.asyn.line.int2,mod.asyn.line.int3)

mod.asyn.line.int <- lmer(asyn_roi_line_avg_int^(1/6) ~ 1 + Groups.f + (1|attempt.no), data = ens4)

# Develop Avg Int a-Syn model designed for back-transformation
warp.ens.asyn.avg.int.line <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_roi_line_avg_int) ~ Groups.f + 
        (1|attempt.no) + (1|attempt.no:well)
         , data = ens4))
# Statistical results summary
summary(em.asyn.line.int <- emmeans(warp.ens.asyn.avg.int.line, ~Groups.f, type = "response"))
pwpm(em.asyn.line.int)

# ENS Neurite PFF Models
mod.pff.line.int0 <- lm(pff_roi_line_avg_int^(1/6) ~ 1 + Groups.f, data = ens4)
mod.pff.line.int1 <- lmer(pff_roi_line_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no), data = ens4)
mod.pff.line.int2 <- lmer(pff_roi_line_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|well), data = ens4)
mod.pff.line.int3 <- lmer(pff_roi_line_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no) + (1|attempt.no:well), data = ens4)

# Identify Best Models
lrtest(mod.pff.line.int0,mod.pff.line.int1,mod.pff.line.int2,mod.pff.line.int3)
anova(mod.pff.line.int1,mod.pff.line.int2,mod.pff.line.int3)

# Develop Avg Int PFF model designed for back-transformation
mod.pff.line.int <- lmer(pff_roi_line_avg_int^(1/6) ~ 1 + Groups.f + (1|attempt.no), data = ens4)
warp.ens.pff.avg.int.line <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_roi_line_avg_int) ~ Groups.f + 
        (1|attempt.no) + (1|attempt.no:well)
         , data = ens4))

# Statistical Summary
summary(em.pff.line.int <- emmeans(warp.ens.pff.avg.int.line, ~Groups.f, type = "response"))
pwpm(em.pff.line.int)

# Plot values
plot.ens.asyn.neurite.int <- ggplot(ens4, aes(x=Groups.f, y=asyn_roi_line_avg_int, fill=Groups.f)) +
  ylab("Average a-Syn Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Neurite a-Syn Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("LPS", "Control")), annotations="**",
              y_position = 50, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("LPS", "But+LPS")), annotations="**",
              y_position = 50, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "6XPFF")), annotations="*",
              y_position = 60, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFFl", "PFF+But")), annotations="*",
              y_position = 50, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = ens4, aes(x = Groups.f, y = asyn_roi_line_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  scale_fill_viridis_d()

plot.ens.asyn.neurite.int

plot.ens.pff.neurite.int <- ggplot(ens4, aes(x=Groups.f, y=pff_roi_line_avg_int, fill=Groups.f)) +
  ylab("Average PFF (ATTO-590) Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Neurite PFF Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 60, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "2XPFF")), annotations="***",
              y_position = 75, tip_length = 0.02, size = 0.8, textsize=8, vjust=0) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="***",
              y_position = 120, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="***",
              y_position = 120, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "6XPFF")), annotations="***",
              y_position = 135, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = ens4, aes(x = Groups.f, y = pff_roi_line_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,145)) +
  scale_fill_viridis_d()

plot.ens.pff.neurite.int

em.ens.asyn.line.int<- plot(em.asyn.line.int, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Neurite a-Syn Average Intensity",
       x = "Estimated Marginal Mean (average intensity a-Syn)",
       y = "Groups")
em.ens.asyn.line.int

em.ens.pff.line.int<- plot(em.pff.line.int, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Neurite a-Syn Average Intensity",
       x = "Estimated Marginal Mean (average intensity PFF)",
       y = "Groups")
em.ens.pff.line.int
```

```{r ens optimization soma}
# ENS Area Soma Non-parametric tests
mod.asyn.p.soma.area <- kruskal.test(asyn_p_area_soma ~ Groups.f, data = ens4)
ph.asyn.p.soma.area <- pairwise.wilcox.test(ens4$asyn_p_area_soma, ens4$Groups.f,p.adjust.method = "BH")
ph.asyn.p.soma.area

mod.pff.p.soma.area  <- kruskal.test(pff_p_area_soma ~ Groups.f, data = ens4)
ph.pff.p.soma.area <- pairwise.wilcox.test(ens4$pff_p_area_soma, ens4$Groups.f,p.adjust.method = "BH")
ph.pff.p.soma.area

# Plot values
plot.ens.asyn.soma.area <- ggplot(ens4, aes(x=Groups.f, y=sqrt(asyn_p_area_soma), fill=Groups.f)) +
  ylab("Sqrt fraction area within soma (a-Syn)") + 
  xlab("Groups") + 
  labs(
    title = "Soma a-Syn Area") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
  geom_signif(comparisons=list(c("Control", "4XPFF")), annotations="**",
              y_position = 0.025, tip_length = 0.0025, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="*",
              y_position = 0.015, tip_length = 0.0025, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "6XPFF")), annotations="***",
              y_position = 0.035, tip_length = 0.0025, size = 0.8, textsize=8, vjust=0.4) +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,0.08)) +
  scale_fill_viridis_d()

plot.ens.asyn.soma.area

plot.ens.pff.soma.area <- ggplot(ens4, aes(x=Groups.f, y=sqrt(pff_p_area_soma), fill=Groups.f)) +
  ylab("Sqrt fraction area within soma (PFF tag)") + 
  xlab("Groups") + 
  labs(
    title = "Soma PFF Area") +  
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 0.27, tip_length = 0.025, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "PFF")), annotations="***",
              y_position = 0.43, tip_length = 0.025, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="***",
              y_position = 0.65, tip_length = 0.025, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+LPS")), annotations="***",
              y_position = 0.27, tip_length = 0.025, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="***",
              y_position = 0.35, tip_length = 0.025, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="***",
              y_position = 0.65, tip_length = 0.025, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("6XPFF", "4XPFF")), annotations="***",
              y_position = 0.75, tip_length = 0.025, size = 0.8, textsize=8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,0.85)) +
  scale_fill_viridis_d()

plot.ens.pff.soma.area

#Avg Int Soma Models
mod.asyn.soma.int0 <- lm(asyn_roi_soma_avg_int^(1/6) ~ 1 + Groups.f, data = ens4)
mod.asyn.soma.int1 <- lmer(asyn_roi_soma_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no), data = ens4)
mod.asyn.soma.int2 <- lmer(asyn_roi_soma_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|well), data = ens4)
mod.asyn.soma.int3 <- lmer(asyn_roi_soma_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no) + (1|attempt.no:well), data = ens4)

# Identify Best Models
lrtest(mod.asyn.soma.int0,mod.asyn.soma.int1,mod.asyn.soma.int2,mod.asyn.soma.int3)
anova(mod.asyn.soma.int1,mod.asyn.soma.int2,mod.asyn.soma.int3)

# Develop Avg Int a-Syn model designed for back-transformation
warp.ens.asyn.avg.int.soma <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_roi_soma_avg_int) ~ Groups.f + 
        (1|attempt.no) + (1|attempt.no:well)
         , data = ens4))

# Statistical results summary
summary(em.asyn.soma.int <- emmeans(warp.ens.asyn.avg.int.soma, ~Groups.f, type = "response"))
pwpm(em.asyn.soma.int)

# PFF ENS Soma Models
mod.pff.soma.int0 <- lm(pff_roi_soma_avg_int^(1/6) ~ 1 + Groups.f, data = ens4)
mod.pff.soma.int1 <- lmer(pff_roi_soma_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no), data = ens4)
mod.pff.soma.int2 <- lmer(pff_roi_soma_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|well), data = ens4)
mod.pff.soma.int3 <- lmer(pff_roi_soma_avg_int^(1/6) ~ 1 + Groups.f + 
                            (1|attempt.no) + (1|attempt.no:well), data = ens4)

# Identify optimal model
lrtest(mod.pff.soma.int0,mod.pff.soma.int1,mod.pff.soma.int2,mod.pff.soma.int3)
anova(mod.pff.soma.int1,mod.pff.soma.int2,mod.pff.soma.int3)

# Develop Avg Int PFF model designed for back-transformation
warp.ens.pff.avg.int.soma <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_roi_soma_avg_int) ~ Groups.f + 
        (1|attempt.no) + (1|attempt.no:well)
         , data = ens4))
summary(em.pff.soma.int <- emmeans(warp.ens.pff.avg.int.soma, ~Groups.f, type = "response"))
pwpm(em.pff.soma.int)

# Plot Values
plot.ens.asyn.soma.int <- ggplot(ens4, aes(x=Groups.f, y=asyn_roi_soma_avg_int, fill=Groups.f)) +
  ylab("Average a-Syn Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Soma a-Syn Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="***",
              y_position = 100, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = ens4, aes(x = Groups.f, y = asyn_roi_soma_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  scale_fill_viridis_d()

plot.ens.asyn.soma.int

plot.ens.pff.soma.int <- ggplot(ens4, aes(x=Groups.f, y=pff_roi_soma_avg_int, fill=Groups.f)) +
  ylab("Average PFF (ATTO-590) Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Soma PFF Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 100, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+LPS")), annotations="*",
              y_position = 100, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "2XPFF")), annotations="***",
              y_position = 115, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="***",
              y_position = 135, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="***",
              y_position = 135, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("6XPFF", "4XPFF")), annotations="***",
              y_position = 150, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="***",
              y_position = 135, tip_length = 0.03, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = ens4, aes(x = Groups.f, y = pff_roi_soma_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,160)) +
  scale_fill_viridis_d()

plot.ens.pff.soma.int

em.ens.asyn.soma.int<- plot(em.asyn.soma.int, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Soma a-Syn Average Intensity",
       x = "Estimated Marginal Mean (average intensity a-Syn)",
       y = "Groups")
em.ens.asyn.soma.int

em.ens.pff.soma.int<- plot(em.pff.soma.int, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Soma a-Syn Average Intensity",
       x = "Estimated Marginal Mean (average intensity PFF)",
       y = "Groups")
em.ens.pff.soma.int
```


```{r ens optimization neuron}
# ENS Area Neuron Non-parametric test
mod.asyn.p.neuron.area <- kruskal.test(asyn_p_area_tot ~ Groups.f, data = ens4)
ph.asyn.p.neuron.area <- pairwise.wilcox.test(ens4$asyn_p_area_tot, ens4$Groups.f,p.adjust.method = "BH")
ph.asyn.p.neuron.area

mod.pff.p.neuron.area  <- kruskal.test(pff_p_area_tot ~ Groups.f, data = ens4)
ph.pff.p.neuron.area <- pairwise.wilcox.test(ens4$pff_p_area_tot, ens4$Groups.f,p.adjust.method = "BH")
ph.pff.p.neuron.area

# Plot Values
plot.ens.asyn.neuron.area <- ggplot(ens4, aes(x=Groups.f, y=sqrt(asyn_p_area_tot), fill=Groups.f)) +
  ylab("Sqrt fraction area within neuron (a-Syn") + 
  xlab("Groups") + 
  labs(
    title = "Neuron a-Syn Area") +
  coord_cartesian(ylim = c(0,0.05)) +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
  geom_signif(comparisons=list(c("Control", "But")), annotations="*",
              y_position = 0.001, tip_length = 0.002, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "4XPFF")), annotations="*",
              y_position = 0.008, tip_length = 0.002, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "6XPFF")), annotations="0.0522",
              y_position = 0.015, tip_length = 0.002, size = 0.8, textsize=4, vjust=0) +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) +
  scale_fill_viridis_d()

plot.ens.asyn.neuron.area

plot.ens.pff.neuron.area <- ggplot(ens4, aes(x=Groups.f, y=sqrt(pff_p_area_tot), fill=Groups.f)) +
  ylab("Sqrt fraction area within neuron (PFF tag)") + 
  xlab("Groups") + 
  labs(
    title = "Neuron PFF Area") +  
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "BLACK") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 0.2, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But")), annotations="0.0581",
              y_position = 0.2, tip_length = 0.02, size = 0.8, textsize=4, vjust=0) +
  geom_signif(comparisons=list(c("PFF", "PFF+LPS")), annotations="**",
              y_position = 0.27, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="***",
              y_position = 0.34, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "PFF")), annotations="**",
              y_position = 0.41, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="***",
              y_position = 0.5, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="***",
              y_position = 0.5, tip_length = 0.02, size = 0.8, textsize=8, vjust=0.4) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,0.6)) +
  scale_fill_viridis_d()

plot.ens.pff.neuron.area

# Develop Avg Int as-Syn model designed for back-transformation
mod.asyn.neuron.int <- lmer(asyn_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + 
                              (1|attempt.no) + (1|attempt.no:well), data = ens4)
warp.ens.asyn.avg.int.neuron <- with(fixed.avg.int.tran, 
    lmer(linkfun(asyn_roi_tot_avg_int) ~ Groups.f + 
        (1|attempt.no) + (1|attempt.no:well)
         , data = ens4))

# Summarize Statistical Results
summary(em.asyn.neuron.int <- emmeans(warp.ens.asyn.avg.int.neuron, ~Groups.f, type = "response"))
pwpm(em.asyn.neuron.int)

# Develop Avg Int PFF model designed for back-transformation
mod.pff.neuron.int <- lmer(pff_roi_tot_avg_int^(1/6) ~ 1 + Groups.f + (1|attempt.no), data = ens4)
warp.ens.pff.avg.int.neuron <- with(fixed.avg.int.tran, 
    lmer(linkfun(pff_roi_tot_avg_int) ~ Groups.f + 
        (1|attempt.no)
         , data = ens4))

# Summarize Statistical Results
summary(em.pff.neuron.int <- emmeans(warp.ens.pff.avg.int.neuron, ~Groups.f, type = "response"))
pwpm(em.pff.neuron.int)

# Plot Values
plot.ens.asyn.neuron.int <- ggplot(ens4, aes(x=Groups.f, y=asyn_roi_tot_avg_int, fill=Groups.f)) +
  ylab("Average a-Syn Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Neuron a-Syn Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("6XPFF", "Control")), annotations="*",
              y_position = 93, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("LPS", "Control")), annotations="*",
              y_position = 85, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("LPS", "But+LPS")), annotations="**",
              y_position = 85, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="*",
              y_position = 85, tip_length = 0.05, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = ens4, aes(x = Groups.f, y = asyn_roi_tot_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) + 
  coord_cartesian(ylim = c(0,100)) +
  scale_fill_viridis_d()

plot.ens.asyn.neuron.int

plot.ens.pff.neuron.int <- ggplot(ens4, aes(x=Groups.f, y=pff_roi_tot_avg_int, fill=Groups.f)) +
  ylab("Average PFF (ATTO-590) Intensity") + 
  xlab("Groups") + 
  labs(
    title = "Neuron PFF Average Intensity") +
  geom_bar(position = 'dodge', stat = 'summary', fun = 'mean') +
    stat_summary(fun.data = mean_se,  
                 geom = "errorbar", position = 'dodge', colour = "black") + 
  geom_signif(comparisons=list(c("PFF", "Control")), annotations="***",
              y_position = 55, tip_length = 0.04, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("Control", "2XPFF")), annotations="***",
              y_position = 70, tip_length = 0.04, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "4XPFF")), annotations="***",
              y_position = 105, tip_length = 0.04, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("2XPFF", "6XPFF")), annotations="***",
              y_position = 120, tip_length = 0.04, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("4XPFF", "4XPFF+But")), annotations="***",
              y_position = 105, tip_length = 0.04, size = 0.8, textsize=8, vjust=0.4) +
  geom_signif(comparisons=list(c("PFF", "PFF+But+LPS")), annotations="***",
              y_position = 55, tip_length = 0.04, size = 0.8, textsize=8, vjust=0.4) +
  geom_dotplot(data = ens4, aes(x = Groups.f, y = pff_roi_tot_avg_int),
               alpha = .8, position_dodge(width = 0.9), binaxis = "y", dotsize = 0.1,
               stackdir = "center", stackratio = 1) +
  theme(legend.position="none",
        text = element_text(size = 14),
        axis.text.x = element_text(angle = 90),
        axis.title.x = element_blank()) +
  coord_cartesian(ylim = c(0,130)) +
  scale_fill_viridis_d()

plot.ens.pff.neuron.int

em.ens.asyn.neuron.int<- plot(em.asyn.neuron.int, 
                                type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Neuron a-Syn Average Intensity",
       x = "Estimated Marginal Mean (average intensity a-Syn)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())

em.ens.asyn.neuron.int

em.ens.pff.neuron.int<- plot(em.pff.neuron.int, type = "response", comparisons = TRUE, plotIT=FALSE) + 
  labs(title = "Neuron a-Syn Average Intensity",
       x = "Estimated Marginal Mean (Average Intensity PFF)",
       y = "Groups") +
  theme(text = element_text(size = 14),
          axis.title.y = element_blank())

em.ens.pff.neuron.int
```


